
###################################
## Read in simulations
###################################
all_sims <- read.csv('data/final_dp_sims.csv', stringsAsFactors = F)
dim(all_sims)

color_vec <- c('orange','navyblue')

# Epsilon plot 
e <- all_sims[which(all_sims$N == 100000 & all_sims$alpha == 0.25 & all_sims$epsilon <= 1 & all_sims$P == 1000),] %>% 
  group_by(epsilon) %>% 
  summarise(tt = mean(theta_tilde),th = mean(theta_hat), n = n())
e <- gather(e, 'type','est',tt:th)

e$bias <- e$est - 3
  
p1 <- ggplot(e) + geom_line(aes(x = epsilon, y = bias, group = type, colour = type)) + 
  geom_point(aes(x = epsilon, y = bias, group = type, colour = type)) + theme_bw() + 
   geom_hline(aes(yintercept = 0), linetype = 'dashed') + 
  labs(x = TeX('$\\epsilon$'), y = 'Bias', colour = '') + 
  scale_colour_manual(labels = c(expression(hat(theta)[dp]), expression(tilde(theta)[dp])), values = color_vec) + 
  annotate("text", x = 0.4, y = -0.008, label = expression(tilde(theta)^{dp}), colour = color_vec[2]) + 
  annotate("text", x = 0.4, y = -0.06, label = expression(hat(theta)^{dp}), colour = color_vec[1]) + 
  guides(colour = FALSE) +
  scale_x_continuous(breaks = c(0.1, 0.25, .5, 0.75, 1))
p1

###################################
# Alpha plot 
###################################
e <- all_sims[which(all_sims$N == 100000 & all_sims$alpha <= 0.75 & all_sims$epsilon == 1 & all_sims$P == 1000),] %>% 
  group_by(alpha) %>% 
  summarise(tt = mean(theta_tilde),
            th = mean(theta_hat),
            n = n())

e <- gather(e, 'type','est',tt:th)
e$bias <- e$est - 3

p2 <- ggplot(e) + geom_line(aes(x = alpha, y = bias, group = type, colour = type)) + 
  geom_point(aes(x = alpha, y = bias, group = type, colour = type)) + theme_bw() + 
  geom_hline(aes(yintercept = 0), linetype = 'dashed') + 
  labs(x = TeX('$\\alpha_2$'), y = 'Bias', colour = '') + 
  scale_colour_manual(labels = c(expression(hat(theta)[dp]), expression(tilde(theta)[dp])), values = color_vec) + 
  scale_x_continuous(labels = unique(e$alpha), breaks = unique(e$alpha)) +
  annotate("text", x = 0.29, y = -0.15, label = expression(hat(theta)^{dp}), colour = color_vec[1]) + 
  annotate("text", x = 0.48, y = -0.03, label = expression(tilde(theta)^{dp}), colour = color_vec[2]) + 
  guides(colour = FALSE)
p2


###################################
# N size plot 
###################################
e <- all_sims %>% filter(alpha == 0.25, epsilon == 1 & ((N != 100000) | (N == 100000 & P == 1000))) %>% 
  group_by(N, alpha, epsilon) %>% 
  summarise(tt = mean(theta_tilde),
            th = mean(theta_hat_dp_before_fix),
            n = n())

e2 <- all_sims[all_sims$alpha == 0.25 & all_sims$epsilon == 0.1 & all_sims$N > 100000,] %>% 
  group_by(N, alpha, epsilon) %>% 
  summarise(tt = mean(theta_tilde),
            th = mean(theta_hat),
            n = n())
e <- rbind(e, e2)

e <- gather(e, 'type','est',tt:th)

e <- e %>% filter(N > 7500)
e$bias <- e$est - 3

p3 <- ggplot(e) + geom_line(aes(x = N, y = bias, group = type, colour = type)) + 
  geom_point(aes(x = N, y = bias, group = type, colour = type)) + theme_bw() + 
  geom_hline(aes(yintercept = 0), linetype = 'dashed') + 
  labs(x = TeX('N'), y = 'Bias', colour = '') + 
  scale_colour_manual(labels = c(expression(hat(theta)[dp]), expression(tilde(theta)[dp])), values = color_vec) +
  annotate("text", x = 500000, y = -0.055, label =  expression(hat(theta)^{dp}), colour = color_vec[1]) + 
  annotate("text", x = 700000, y = -0.015, label = expression(tilde(theta)^{dp}), colour = color_vec[2]) + 
  scale_x_continuous(breaks = c(10000, 250000, 500000, 750000, 1000000))+
  guides(colour = FALSE) + theme(axis.text.x = element_text(size = 7))
p3



#########################################
# Variance plot with epsilon  0.1 to 1
#########################################
e <- all_sims[which(all_sims$N == 100000 & all_sims$alpha == 0.25 & all_sims$epsilon <= 1 & all_sims$P == 1000),] %>% 
  group_by(epsilon) %>% 
  summarise(est_sd = sqrt(mean(var_est)),
            sim_sd = sqrt(var(theta_tilde)),
            sd_theta_hat = sqrt(var(theta_hat)),
            n = n())
e <- gather(e, 'type','est', est_sd:sd_theta_hat)



p4 <- ggplot(e) + geom_line(aes(x = epsilon, y = est, group = type, colour = type)) + 
  geom_point(aes(x = epsilon, y = est, group = type, colour = type)) + theme_bw() + 
  geom_hline(aes(yintercept = 0), linetype = 'dashed') + 
  labs(x = TeX('$\\epsilon$'), y = 'Std. Error', colour = '') + 
  scale_colour_manual( values = c('navyblue','orange', 'grey40')) +
  annotate("text", x = 0.4, y = 0.09, label = expression(paste(SE[hat(theta)^dp])), colour = 'orange') +
  annotate("text", x = 0.2, y = .03, label = expression(paste(hat(SE)[tilde(theta)^{dp}])), colour = 'navyblue') + 
  annotate("text", x = 0.75, y = 0.015, label = expression(paste(SE[tilde(theta)^{dp}])), colour = 'grey40') + 
  guides(colour = FALSE) + 
  scale_x_continuous(breaks = c(0.1, 0.25, .5, 0.75, 1))
p4


###################################
# P plot 
###################################
e <- all_sims %>% 
  filter(epsilon == 1, alpha == 0.25, N == 100000) %>% 
  #filter(alpha == 0.25, N == 100000) %>% 
  group_by(P) %>% 
  summarise(tt = mean(theta_tilde),
            th = mean(theta_hat_dp_before_fix),
            n = n(),
            N = mean(N),
            epsilon = mean(epsilon))
e <- e%>% filter(P <= 500)
e <- gather(e, 'type','est', tt:th)
e$bias <- e$est - 3

p5 <- ggplot(e) + geom_line(aes(x = P, y = bias, group = type, colour = type)) + 
  geom_point(aes(x = P, y = bias, group = type, colour = type)) + theme_bw() + 
  geom_hline(aes(yintercept = 0), linetype = 'dashed') + 
  labs(x = TeX('P'), y = 'Bias', colour = '') + 
  scale_colour_manual(labels = c(expression(hat(theta)[dp]), expression(tilde(theta)[dp])), values = color_vec) +
  annotate("text", x = 50, y = 0.2, label =  expression(hat(theta)^{dp}), colour = color_vec[1]) + 
  annotate("text", x = 250, y = -0.03, label = expression(tilde(theta)^{dp}), colour = color_vec[2]) + 
  #scale_x_continuous(breaks = c(200, 1000, 2000, 3000, 4000))+
  guides(colour = FALSE)
p5


###################################
# Save plots
###################################

bigplot <- ggarrange(p2, p1, p3, p4)
bigplot
ggsave(bigplot, file = 'plots/fig3.pdf', width = 7, height = 5)

ggsave(p5, file = 'plots/fig4.pdf', width = 7, height = 5)
